How input fluctuations reshape the dynamics of a biological switching system 
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An important task in quantitative biology is to understand the role of stochasticity in biochemical 
regulation. Here, as an extension of our recent work [Phys. Rev. Lett. 107, 148101 (2011)], we 
study how input fluctuations affect the stochastic dynamics of a simple biological switch. In our 
model, the on transition rate of the switch is directly regulated by a noisy input signal, which 
is described as a nonnegative mean-reverting diffusion process. This continuous process can be a 
good approximation of the discrete birth-death process and is much more analytically tractable. 
Within this new setup, we apply the Feynman-Kac theorem to investigate the statistical features of 
the output switching dynamics. Consistent with our previous findings, the input noise is found to 
effectively suppress the input-dependent transitions. We show analytically that this effect becomes 
significant when the input signal fluctuates greatly in amplitude and reverts slowly to its mean. 

PACS numbers: 02.50.Le, 05.65. +b, 87.23.Ge, 87.23.Kg 



I. INTRODUCTION 

Stochasticity appears to be a hallmark of many bio- 
logical processes involved in signal transduction and gene 
regulation [TlfTTj . Over the past decade, there have been 
numerous experimental and theoretical efforts to under- 
stand the functional roles of noise in various living sys- 
tems 2 41 . Remarkably, the building block of different 
regulatory programs is often a simple two-state switch 
under the regulation of some noisy input signal. For ex- 
ample, a gene network is composed of many interacting 
genes, each of which is a single switch regulated by spe- 
cific transcription factors. At a synapse, switching of 
ligand-gated ion channels are responsible for converting 
the presynaptic chemical message into postsynaptic elec- 
trical signals, while the opening and closing of each chan- 
nel depends on the binding of certain ligands (such as a 
neurotransmitter). In bacterial chemotaxis, the cellular 
motion is powered by multiple flagellar motors and each 
motor rotates clockwise or counterclockwise depending 
on the level of specific response regulator (e.g., CheY-P 
in E. coli). Irrespective of the context, the input sig- 
nal, i.e., the number of regulatory molecules, is usually 
stochastic due to discreteness, diffusion, random birth 
and death. How does a biological switching system work 
in a noisy environment? This is the central question we 
attempt to address in this paper. 

Previous studies on this topic have usually focused on 
the approximate, static relationship between input and 
output variations (e.g., the additive noise rule) [T2lU9| . 
while the dynamic details (e.g., dwell time statistics) of 
the switching system have often been ignored. Our re- 
cent work suggests that there is more to comprehend 
even in the simplest switching system |22j . For example, 
we showed that increasing input noise does not always 
lead to an increase in the output variation, disagreeing 
with the additive noise rule as derived from the coarse- 



grained Langevin approach. Traditional methods often 
use a single Langevin equation to approximate the joint 
input-output process, and relies on the assumption that 
the input noise is small enough such that one can lin- 
earize input-dependent nonlinear reaction rates. Our ap- 
proach to this problem is quite different as we explicitly 
model how the input stochastic process drives the output 
switch, without making any small noise assumption. 

In our previous paper (55], the input signal was gen- 
erated from a discrete birth-death process and regulated 
the on transition of a downstream switch. By explicitly 
solving the joint master equation of the system, we found 
that input fluctuations can effectively reduce the on rate 
of the switch. In this paper, this problem is revisited 
in a continuous noise formulation. We propose to model 
the input signal as a general diffusion process, which is 
mean- reverting, nonnegative, and tunable in Fano factor 
(the ratio of variance to mean). We employ the Feynman- 
Kac theorem to calculate the input-dependent dwell time 
distribution and examine its asymptotic behavior in dif- 
ferent scenarios. Within this new framework, we recover 
several findings reported in [32], and also demonstrate 
how the noise-induced suppression depends on the rel- 
ative noise level as well as the relative input-to-output 
timescale. Finally, we elaborate on how the diffusion 
process introduced in this paper can be a reasonable ap- 
proximation of the discrete birth-death process. 



II. MODEL 

The input of our model, denoted by X(t), represents 
a specific chemical concentration at time t and directly 
governs the transition rates of a downstream switch. The 
binary on-off states of the switch in continuous time con- 
stitute the output process, Y(t). A popular choice for 
X(t) is the Ornstein-Uhlenbeck (OU) process due to its 
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FIG. 1: Illustration of our model: X(t) represents the input 
signal which fluctuates around a mean level over time; Y(t) 
records the switch process which flips between the off (Y = 0) 
and on (Y = f) states with transition rates k on X(t) and fc ff • 
t is the dwell time in the off state. 



analytical simplicity and mean-reverting property [^2T - 
|4"4"] . However, this process does not rule out negative 
values, an unphysical feature for modeling chemical con- 
centrations. For both mathematical convenience and bio- 
physical constraints (see Section IV for more details), we 
model X(t) by a square-root diffusion process 



dX(t) = X[fi - X{t)}dt + ayJX{t)dW u 



(1) 



where A represents the rate at which X(t) reverts to its 
mean level fi, a controls the noise intensity, and Wt de- 
notes the standard Brownian motion. This simple pro- 
cess is known as the Cox-Ingersoll-Ross (CIR) model of 
interest rates [35]. The square-root noise term not only 
ensures that the process never goes negative but also cap- 
tures a common statistical feature underlying many bio- 
chemical processes, that is, the standard deviation of the 
copy number of molecules scales as the square root of 
the copy number, as dictated by the Central Limit Theo- 
rem. Solving the Fokker-Planck equation for Eq. (1), we 
obtain the steady-state distribution of the input signal: 
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one can also find the steady-state covariance: 



Cov[X(t),X(t + s)\ 
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Thus X(t) is a 



stationary process with correlation time A -1 . 

In reality, the output switching rates may depend on 
the input X(t) in complicated ways, depending on the de- 
tailed molecular mechanism. For analytical convenience, 
we assume that the on and off transition rates of the 
switch are k on X(t) and k g, respectively (Fig. I). As 



a result, the input fluctuations should only affect the 
chance of the switch exiting the off state; the on-state 
dwell time distribution is always exponential with rate 
parameter k Q g. If we mute the input noise, then Y(t) re- 
duces to a two-state Markov process with transition rates 
k on fi and k g. However, the presence of input noise will 
generally make Y(t) a non-Markovian process, because 
the off-state time intervals may exhibit a non-exponential 
distribution. To illustrate this point more rigorously, we 
consider the following first passage time problem [44] . 
Suppose the switch starts in the off state at t = with 
the initial input X(0) = x. Let r be the first time of the 
switch turning on (Fig. 1). Then the survival probability 
f(x, t) for the switch staying off up to time t is given by 



f(x, t) = P(t > t\X(0) = x) = E 1 



- /* k oa x{s)ds\ i (3) 



where E a [...] denotes expectation over all possible sample 
paths of X(s) for < s < t, conditioned on X(0) = x. 
The Feynman-Kac formula [44] asserts that f(x,t) must 
solve the following partial differential equation: 
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(4) 



with the initial condition f(x,0) — 1. As we will show 
in the next section, the off-state dwell time distribution 
is not exactly exponential, though it is asymptotically 
exponential when t is large. 



III. RESULTS 

A similar partial differential equation to Eq. Q has 
been solved in Ref. [45j to price zero-coupon bonds under 
the CIR interest rate model. The closed-form solution for 
our problem is similar and is found to be: 



f(x,t) 
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which is a Gamma distribution with stationary variance where 
o\ = /icr 2 /(2A). This is another attractive aspect of 
this model, since the protein abundance from gene ex- 
pression experiments can often be fitted by a Gamma 
distribution [8HTU]. The parameter a in Eq. (2) can be 
interpreted as the signal-to-noise ratio, since a = /i 2 jo\. 
For a > 1, the zero point is guaranteed to be inaccessible 
for X(t). The other shape parameter f3 sets the Fano 
factor as we have 1//3 = o\j[L. Using the Ito calcu- 
lus 
linit 



A = VA 2 + 2fc on( 7 2 = Av/1 + 2fc on a 2 /A 2 



(5) 
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Evidently, f(x,t) remembers the initial input x and de- 
cays with t in a manner which is not exactly exponential. 
For t A -1 , Eq. Q5J) takes the following form, 
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In deriving Eq. Q, we have used the following relation- 
ship which defines the new parameter k on as below: 



= q(A — A) 
= 2/x 



A 



(A-A) = 



2A 
A + A 



(8) 



Thus, asymptotically speaking, f(x,t) decays with t in 
an exponential manner at the rate k on fj,. It can be easily 
verified that Eq. ^ is a particular solution to Eq. Q . 

To gain further insight into how the input noise affects 
the switching dynamics, we first study the "slow switch" 
limit where X(t) fluctuates so rapidly that A -1 <C Ty- 
Here, Ty = (fconM + fcoff) -1 is the output correlation time 
for the noiseless input model. In this limit, the initial 
input x in f(x,t) at the start of each off period is effec- 
tively drawn from the Gamma distribution P s (x) defined 
in Eq. (2); the successive off-state intervals are almost 
independent with each other (as the input loses memory 
quickly) and are distributed as 



P(r < t) = 1 - P(t > t) = 1 



f(x,t)P s (x)dx. (9) 



By direct integration over x, we find that 
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exp(-fc on ^). (10) 



In the last step above, we have used Eq. (8) as well as 
the following observation: By introducing 9 = k on a 2 /A 2 
which reflects the deviation of A from A, one can check 
that as long as A ^ fc OI1 /j (as ensured by A -1 <C Ty) the 
following holds, regardless of the values of 9, 
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Our simulations show that the approximate result, P(r > 
t) ~ g-fcon^t^ j g excellent (Fig. 2A), independent of the 
values of 9. Thus, the average waiting time for the switch 
to turn on is (fconA*) 1 ; longer than the corresponding 
average time (fc ony it) -1 for the noiseless input model. This 
is similar to our previous result and suggests that 
the input noise will effectively suppress the on state by 
increasing the average waiting time to exit the off state. 
Consequently, the probability, P on , to find the switch on 
(Y = 1) is less than that in the noiseless input model: 
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where K d = k s/k on , the effective equilibrium constant, 
is larger than the original K d = k g/k on as per Eq. (8). 

Therefore, in the slow switch limit (Ty ^S> A -1 ), the 
output Y(t) is approximately a two-state Markov process 
with transition rates k on ii and fc ff- If we further assume 
that the noise is modest (i.e., ax < ff), then 
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FIG. 2: (color online). The slow switch case. Here we use 
A = 10, fcon = 0.02, and k oS = 0.1 (thus K d = 5). (A) P{t > 
t) versus t for /i = 3 and 8 = 0.005, 0.5, and 2.0 which are 
achieved by choosing a = 5, 50, and 100. Symbols represent 
simulation results, while lines denote exp(— k on /j,t). (B) a Y — 
a Y versus ax, where different values of a x are obtained by 
tuning a. (C) a Y and a\ versus [i/Kd with 8 = 0.50. (D) 
Ty — Ty versus cj\ . 



The equality above shows that 9 is a characteristic pa- 
rameter determined by the ratio of the input to output 
time scales and the relative noise strength. For 9 <C 1, 
we have A = Ayl + 26 ~ A(l + 9) by Eq. (6), and the 
effective on rate defined in Eq. (8) becomes, 
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In [35] , the input signal was taken to be a Poisson birth- 
death process for which the variance is equal to the mean 
(<j x = ff), and the time scale has been normalized by 
putting the death rate equal to one (which amounts to 
setting A = 1 here). These two constraints reduce Eq. 
(13) to k on ~ fc on /(l + fcon), recovering the result we 
derived in 22 . The consistency indicates that our key 
findings are general, independent of the specific model we 
choose. The continuous diffusion model here, however, is 
more flexible as it allows the Fano factor to differ from 
one, that is, a\ ^ /i. 

For small 9, the stationary variance of Y(t) can be 
expanded as follows: 



= P on (l - P on ) ~ a 2 Y 
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FIG. 3: (color online). The fast switch case. Here we choose 
A = 0.0125, fi = 5, k oS = 0.1, fc on = 0.02, and 9 = 10 (such 
that a ~ 0.28). (A) Sample ACF of successive off-state time 
intervals. (B) Distribution of the off-state intervals, P(r > t). 
Symbols are from Monte-Carlo simulations and solid line is 
the semi-analytical prediction described in the main text. (C) 
Sample ACF of a simulated sample path of Y(t) in the semi- 
log scale. (D) Input distributions conditioned on the output. 



where o~ Y = 



fj,Kd/(fJ, + Kd) 2 is the output variance of the 
noiseless input model. Eq. ( 14 ) indicates that the input 
noise a\ does not always contribute positively to the 
output variance a Y - In fact, the contribution is negligible 
when fi is near Kd and even negative for /i < Kd (Fig. 
2B). The explanation is, as we argued before [22 , that 
a two-state switch at any moment is a Bernoulli random 
variable whose variance is always bounded by one quarter 
(Fig. 2C). Finally, with the effective on rate fc on , the 

correlation time of Y(t) becomes Ty = (^onM + ^off) 
and can likewise be expanded as follows: 
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Thus, Ty weakly increases with the input noise in this 
small noise limit (Fig. 2D). 

We now examine the "fast switch" limit where the 
switch flips much faster than the input reverts to its 
mean (Ty < ^'). In this scenario, the initial input 
values {xi,i = 1,2,...} for successive first-passage time 
intervals {Ti,i = 1,2, ...} are correlated due to the slow 
relaxation of X(t). This memory makes the sequence 
{Ti,i — 1,2,...} correlated as well, as confirmed by our 
Monte-Carlo simulations (Fig. 3A). For the same rea- 
son, the autocorrelation function (ACF) of the output 
Y(t) exhibits two exponential regimes (Fig. 3C): over 
short time scales, it is dominated by the intrinsic time 



Ty of the switch; over long time scales, however, it de- 
cays exponentially at the input relaxation rate A. This 
demonstrates that the long-term memory in X(t) is in- 
herited by the output process Y(t). For a fast switch, 
the distribution P(r > t) is not fully exponential (Fig. 
3B), though it decays exponentially at the rate of k on /i 
for large t, as predicted by the asymptotic Eq. (7). We 
can still use the closed-form solution of f(x, t) in Eq. (5) 
to fit the simulation data (open circles). Specifically, we 
calculate Pir > t) — f{x 1 t)Pf(x)dx (blue line in Fig. 
3B), where P?(x) is the distribution of the initial x for 
each switching event (defining the "first-passage" time r) 
and can be obtained from the same Monte-Carlo simu- 
lations. Clearly, such a semi- analytical trial (blue line) 
provides a nice fit to the simulation results (open circles) . 

Note that (x) ^ P s (x) due to memory effects in the 
fast switch limit. To show this, we illustrate the input- 
output interdependence by plotting the input distribu- 
tions conditional on the output state in Fig. 3D. In fact, 
we shall have P(X = x\Y = 1) = Pt(x), because: first, x 
in Pt(x) denotes x = A"(i,j") where to is the last time of 
off-transition; second, AT(t,|) = X(t^) due to continuity 
and Y(tg) = 1 by definition; third, all the on-state in- 
tervals are memoryless with the same Poisson rate k g. 
This simple relation P(X — x\Y = 1) = Pr(x) has been 
confirmed by our simulations (results not shown). It is 
also obvious from Fig. 3D that the expectation value of 
the input conditioned on Y = 1 is larger than that given 
Y = 0. Thus the mean of X(t) should lie in between, 
i.e., E[X\Y = 1] > E[X] > E[X\Y = 0], an interesting 
feature of the fast switch limit [25]. Since f(x,t) is a 
decreasing function of x and the initial input x is likely 
to be larger under P?{x) = P(X = x\Y = 1) than under 
the measure P S {X — x), we should have 



P{r>t) = f 
Jo 



f(x,t)Pr(x)dx 



< / /(x,<)F s ^)dx~e^ fe ^*. (16) 
Jo 

The above inequality explains why the distribution of r is 
below the single exponential e _fcon ' i * (dashed line) in Fig. 
3B. All the above results (Fig. 3A-D) indicate that the 
output process Y(t) is non-Markovian in the fast switch 
limit and, again, confirm the more general applicability 
of our findings reported in [55] . 



IV. DIFFUSION APPROXIMATION 

In this paper, we have used the square-root diffu- 
sion (or CIR) process to model biochemical fluctua- 
tions. Here we argue that this choice is inspired by 
the fundamental nature of general biochemical processes. 
Many biochemical signals are subject to counteracting 
effects: synthesis/degradation, activation/deactivation, 
transport in/out of a cellular compartment, etc. As a 
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result, these signals tend to fluctuate around their equi- 
librium values. A simple yet realistic model to capture 
these phenomena is the birth-death process, which we 
adopted to model biochemical noise in [53] . Remarkably, 
a birth-death process can be approximated by a Markov 
diffusion process [42] . The standard procedure is to em- 
ploy the Kramers-Moyal expansion to convert the master 
equation into a Fokker-Planck equation (if terminating 
after the second term). This connection allows the use 
of a Langevin equation to approximate the birth-death 
process. We will explain this in a more intuitive way. 

Assume that the birth and death rates for the input 
signal X(t) are v and XX (t). Then the stationary distri- 
bution of X(t) is a Poisson distribution, with its mean, 
variance, and skewness given by /i = v/X, o\ = fx, and 
respectively. The corresponding Langevin equa- 



"shift" in distribution. This becomes a limitation of the 



Langevin approximation Eq. (17) which could fail if the 



noise is large (i.e. the number of molecules is small) . 

Under Eq. ( 19 ) for the input X(t), we can still evaluate 
the analog of Eq. (3): 



f(x,t) = W 



- f ( t ) k aa X(s)ds 



-1/2 



= E x+fJ " e~fo k „X'(s)ds e k oat it^ ^21) 

where E^+^f...] denotes expectation over all possible sam- 
ple paths of X' over [0, f], conditioned on X'(0) = x + \l. 
Since X' follows the CIR process Eq. ( p0| , the expec- 
tation E^+^f...] has an expression similar to Eq. (5). In 
fact, the survival probability f(x,t) in Eq. |2l]) equals 



tion that approximates this birth-death process is 



dX(t) 
dt 



v-xx{t) + <q(t), 



(17) 



Ae A '/ 2 /sinh(Ai/2) 
A + Acoth(Af/2) 



exp 
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X + Acoth(Af/2) 



where the stochastic term 77(f) represents a white noise 
with (77(f)) = and delta-correlation 

(?7(f)77(f')) = (v + XX)5{t-t') = X((i + X)S(t-t'). (18) 



with A = VA 2 + 2fc on A and a = 4/i. At large f, this is 



f(x,t) ~ exp 
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(2^on k on )flt 



(22) 



Physically, the Langevin approximation Eq. ([17]) holds 
when the copy number of molecules is large and the time 
scale of interest is longer than the c har acteristic time 
(A -1 ) of the birth-death process. Eq. (18) indicates that 
the instantaneous variance of the noise term 77(f) equals 
the sum of the birth rate v and the death rate XX(t). 
An intuitive interpretation is that since both the birth 
and death events follow independent Poisson processes, 
the total variance of the increment X(t + At) — X(t) over 
a short time Af must be equal to z/Af + XX(t)At. As 
fi = v/X, we can rewrite the Langevin Eq. (17) as the 



where k' on = 2fc on /(l + y 7 ! + 26") and 9' = k on /X. Thus, 
the dwell time distribution behaves asymptotically as 



P{t > t) ~ e- fe °»^, where k oa = 2k' oa - k oa . (23) 
For 9' <§; 1, the asymptotic rate above becomes: 
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The first equality above shows that fc on can be negative 
when 9' > 4 or k on > 4A. This arises as a curse of the 



following Ito-type stochastic differential equation (SDE): possibility that X(t) can go negative under Eq. ([19]) 



dX(t) = X\jx - X(t)]dt + VA/i + XX(t)dW t , (19) 

which is similar to Eq. (1) introduced at the beginning. 
A transformation X'(t) = X(t) + fj, is convenient for ex- 
ploiting our existing results, as the SDE for X'(t) is 



When dealing with the Langevin approximation Eq. 



(17), researchers usually assume that the input noise is 



small enough (given a large copy number) such that the 
random variable X could be replaced by its mean ji in 
Eq. (18). This results in an OU approximation: 



dX' = A(2/x - X')dt + VxX'dWt. 



(20) 



dXit) = A[/i - X{t)]dt + V2vdW t 



(25) 



which is a particular CIR process. It is easy to check 
that the stationary variance of X'(t) equals /i, and so 
does the variance of X(t). In other words, we still have 



which takes a Gaussian distribution in steady state with 
a\ = [i. Compared to the shifted Gamma distribution 
resulted from Eq. (19), the Gaussian model of X(t) has 



o~x = A* for X(t) under Eq. ( 19 ). As a CIR process, X' in 



equilibrium follows a Gamma distribution, the skewness 
of which is found to be /i~ 1//2 . Therefore, X = X' — fx 
follows a "shifted" Gamma distribution with its mean, 
variance, and skewness given by fj,, fi, and /i -1 ^ 2 , respec- 
tively, the same to those of the Poisson distribution. This 
matching of moments suggests that Eq. ( 19 ) is indeed a 



zero skewness and is unbounded below. Thus, when /j is 
small, the OU approximation becomes an inappropriate 
choice. By Feynman-Kac theorem, the survival probabil- 
ity f(x,t) under Eq. (25) satisfies 
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nice approximation to the original birth-death process. 
However, X(t) in Eq. (19) can take negative values, be- 
cause X = X' — \i is bounded below by —\i due to its 



which can also be exactly solved. We omit the solution 
here, but later will show that P(X(t) < 0) > for the 
OU process will lead to k oa < in certain regimes. 
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constant. Inspired by Eqs. ^ and (22), we guess that 
when t is sufficiently large, f(x,t) ~ exp(— Ax — fc on /ii), 
for some constant coefficients A and k on (to be deter- 
mined). Plugging this expression into Eq. ( |29[ ) yields: 

-k on fif = -AX(^-x)f+^A 2 (a 2 ) +a 2 1 x)f-k on xf, (30) 

which holds only when A and fc on jointly solve the fol- 
lowing two algebraic equations: 
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FIG. 4: (color online). Comparison of the Poisson, "shifted" 
Gamma, Gaussian, and Gamma distributions, all of which 
have the same mean fi = 3 and the same variance a x = [l- 



We propose an alternative fix for the Langevin approx- 
imation, replacing the mean [i by its random counterpart 
X(t) in Eq. (pi). This yields a CIR process: 



Given <j\ > 0, the solution of Eq. (32) is: 



A = 



-X + y/X 2 + 2k on a 2 



2k n 



X + yJX 2 + 2k on a 2 ' 



(33) 



Thus, by defining 6' = k oa a 2 /X , Eq. (33) becomes 

9k 



AX 



1 + VI + 26" 



kL 



(34) 



dX(t) = X\p - X(t)]dt + y / 2XX(t)dW t , (27) Eliminating A 2 in Eqs. (31) and (32), we find that 



under which X(t) follows a Gamma distribution in steady 
state, with a x = fi as in all the previous diffusion ap- 
proximations. The skewness in this model is found to 
be 2/z -1 / 2 , which is twice the skewness in the (birth- 
death) Poisson distribution. Though larger than desired, 
it is better than none (in the OU case). Fig. 4 plots 
a comparison of the Poisson, Gamma, shifted Gamma, 
and Gaussian distributions, all satisfying a\ = fi. One 
can see that the shifted Gamma distribution, which fol- 
lows from Eq. (19), gives the closest approximation to 



the (birth-death) Poisson distribution, while the Gamma 
and Gaussian distributions are, roughly speaking, equally 
good to approximate the Poisson. Nonetheless, only the 
Gamma density from the CIR model is nonnegative, like 
the original Poisson distribution. 

The diffusion approximations we have discussed so far, 
including Eq. 0, Eq. @, Eq. J25j), and Eq. ((27}, are 



all special cases of the following general Ito SDE: 



dX(t) = X[n - X{t)]dt + J a 2 + o 2 X(t)dW t . (28) 



Under Eq. (28) the survival probability f{x,t) satisfies 



df 

dt 



1 



,#7 

dx 2 



KnXf. (29) 



Again, a shortcut for solving Eq. (29) is to introduce 
X' = X + (i 2 /a 2 which will evolve as a CIR process. This 
will allow us to make use of the existing results and obtain 
a similar expression for f(x,t) as before. However, our 
main interest is the effective rate fc on in the asymptotic 
behavior of P(r > £) ~ fl exp(— k on )it), where £1 is some 



kn 



-I 



fen 



fiaf 



(35) 



Clearly, when ctq = 0, Eq. (28) becomes a CIR process 
and Eq. S reduces to fc on = k' on = 2fc on /(l + Vl + 29'), 



reco vering our result in Section III. When o 2 /<j 2 = fx, Eq. 



(35) is fc on = 2k' on — fc on , coinciding with Eq. (23) 



Finally, if <j\ — 0, the diffusion Eq. (28) becomes an 
OU process and A = k on /X by Eq. (32). In this case, 



k rr 2 

2^X 2 



fen 



W< ( v 
X 



'x 
H 2 



, (36) 



where a\ = a^/(2X) is the variance of the OU process. 
We can define in the same way as in Eq. ( 12 1, such that 
fcon = k on {l — 8/2) in Eq. (36). This becomes negative 



if 9 > 2, again a consequence of P(X(t) < 0) > 0. Note 
that the OU approximation Eq. (25) is a particular OU 

= fcon(l fcon/ A) 



process with a x — /i. This leads to fc, 
in Eq. |36| 
&on/(l + k on /X), consistent with our result in 



So for small or A > fc on , we have fc on 



In sum, the general diffusion process defined by Eq. 



( 28 1 may take negative values with positive probability, 



unless Oq = which corresponds to the CIR process. 
Negative biochemical input is unrealistic and may lead 
to unphysical results (such as fc on < 0) if the input noise 
is large. For this reason, we conclude that the CIR pro- 
cess Eq. is more suitable for modeling biochemical 
noise in the continuous setup. As shown, it is analyti- 
cally tractable and possesses desirable statistical features, 
including stationarity, mean-reversion, Gamma distribu- 
tion, and a tunable Fano factor. 
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V. CONCLUSION 

In this paper, we have extended our previous research 
on the role of noise in biological switching systems. We 
propose that a square-root diffusion process can be a 
more reasonable model for biochemical fluctuations than 
the commonly used OU process. We employ standard 
tools in stochastic processes to solve a well-defined funda- 
mental biophysical problem. Consistent with our earlier 
results, we find that the input noise acts to suppress the 
input-dependent transitions of the switch. Our analyt- 
ical results in this paper indicate that this suppression 



increases with the input noise level as well as the in- 
put correlation time. The statistical features uncovered 
in this basic problem can provide us with new insights 
to understand various experimental observations in gene 
regulation and signal transduction systems. The current 
modeling framework may also be generalized to incorpo- 
rate other biological features such as ultrasensitivity and 
feedbacks. Work along these lines is underway. 
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